function dy = f(t,y)
    dy = zeros(3,1);
    dy(1) = y(2);
    dy(2) = 1.5 * (1-y(1)^2)*y(2) - y(1);
    dy(3) = y(1)*y(2);
end

